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SUMMARY 


The  primary  results  of  the  Phase  II  efforts  are  the  following: 

1.  Continued  development  and  refinement  of  an  approach  for  real  time 
on-board  estimation  of  spacecraft  attitude  with  sub-five  arc-second 
precision. 

2.  Implementation  and  validation  of  several  variations  of  the  approach 
in  a  laboratory  microcomputer  -  the  objective  being  to  ascess  the 
problems  associated  with  a  real-time,  on-board  version  of  this 
system. 

3.  Development  of  truth  models  to  generate  realistic  input  data  for  the 
star  pattern  recognition  and  Kalman  filter  strategies. 

4.  Conversion  from  use  of  Euler  angles  to  Rodrigues  parameters  to  define 
vehicle  attitude,  affecting  the  algorithms  for  star-pattern  recognition 
least-squares  differential  correction  to  refine  estimated  attitude, 
and  the  Kalman  filter  strategy  to  obtain  the  optimal  attitude  estimate. 

5.  Formulation  of  algorithms  using  Euler  parameters  to  define  orientation. 

These  results  are  discussed  in  detail  herein. 
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1.0  INTRODUCTION  AND  SYSTEM  OVERVIEW 

The  purpose  of  Che  research  activities  described  in  this  document 
is  to  research  and  develop  a  system  for  on-board  attitude  determination 
for  satellites,  achieving  sub-five  arc-second  accuracy  when  such  a 
system  is  interfaced  with  a  general  purpose  star  tracker  (typical  of 
several  existing  configurations). 

The  primary  motivation  for  this  research  is  the  exploitation  of 
recently  developed  light  sensitive  Charge-Doupled-Devices  (CCD)  in  the 
focal  plane  of  a  tracker  lens.  The  CCD  will  act  as  a  "film"  for  imaging 
starlight  and  from  which  the  attitude  can  be  determined.  The  attitude 
determination  can  be  performed  on-board  with  a  high  degree  of  accuracy 
and  high  speed.  A  typical  CCD  for  this  application  consists  of  about 
200,000  light-sensitive  elements  or  pixels,  accurately  imbedded  (to  1 
part  in  10,000)  in  a  microcircuit  chip.  The  CCD  can  be  scanned  electronically 
and  the  voltage  response  of  illuminated  pixels  digitized  for  processing  by 
a  microcomputer;  corrections  can  be  applied  for  image  distortion  based  upon 
apriori  calibration. 

The  image  coordinate  data,  output  by  the  computer,  can  be  either  tele¬ 
metered  to  ground  for  later  analysis  or,  as  described  in  this  report, 
analyzed  on-board  (via  computers  configured  in  parallel)  to  determine 
satellite  attitude  autonomously  in  near  real-time. 

The  basic  system  we  propose  consists  of  2  or  3  CCD  star-trackers 
and  3  microcomputers,  each  with  a  dedicated  function.  The  function  of 
each  of  the  4  sub-systems  is  outlined  below,  with  reference  to  Figures 
1.1,  1.2  and  1.3. 
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1)  CCD  State  Sensors  and  Associated  Electronics. 

Although  the  development  of  CCD  sensors  and  CCD  trackers 
is  not  part  of  this  research,  there  are  several  CCD  star 
tracker  designs  proposed  by  various  organizations  involved 
in  hardware  development.  The  purpose  of  our  work  has  been 
exploitation  of  the  CCD  star  tracker  technology;  we  have 
chosen  a  particular  set  of  parameters  but  it  should  be  kept 
in  mind  that  these  are  nominal,  achievable  values  without 
further  CCD/star  tracker  technology  advances.  Each  of  the  two 
(or  three)  trackers  (assumed  in  various  configurations  herein) 
is  identical  and  their  boresights  are  assumed  separated  by  the 
nominal  interlock  angle  of  90°.  We  assume  a  lens  focal  length 
of  70  mm  and  a  CCD  array  size  such  that  a  7°  x  9°  field  of 
view  (FOV)  is  imaged  onto  the  array.  The  arrays  are  assumed 
to  be  Fairchild's  11.4  mm  x  8.8  mm  matrix  consisting  of  488  x 
380  silicon  pixels.  Starlight  is  defocused  slightly  on  the 
CCD  in  order  to  spread  typical  images  over  9  to  25  adjacent 
pixels.  This  permits  accurate  "centroiding"  of  the  image  to 
determine  image  coordinates  accurate  to  about  10%  of  a  pixel. 
The  processing  of  a  data  frame  consists  of  a  rapid  analog 
readout  of  all  pixels  and  an  analog  to  digital  (A/D)  conver¬ 
sion  only  of  selected  pixels  (based  upon  response  above  an 
analog  threshold  level).  The  scans  of  each  field  of  view  are 
controlled  by  a  common  clock  and  are  assumed  to  represent  2 
(or  3)  frames  (1  from  each  sensor)  taken  at  the  same  instant. 
This  assumption  is  valid  for  all  but  very  rapidly  spinning 
satellites,  since  the  CCDs  can  be  scanned  10  times  per  second. 


2)  Microcomputer  A 

Program  P40C£44  A  is  performed  by  Microcomputer  A  with 
either  one  computer  per  sensor  or  sequential  treatment  of  data 
for  2  or  3  sensors.  Again,  Process  A  is  not  part  of  this 
research  program.  Since  the  functions  to  be  performed  are 
straightforward,  we  simply  replace  Process  A  by  calculating 
synthetic  output  data  whose  availiability  is  clock  controlled. 
Process  A  takes  as  input  data  the  digitized  pixel  voltages 
and  pixel  coordinates  for  up  to  6  stars  in  each  FOV  and  the 
associated  time.  Image  centroids  are  calculated  for  each 
image  and  corrections  for  lens  distortion  and  other  known 
error  sources  are  applied  and  a  relative  magnitude  or  intensity 
is  calculated.  The  output  of  Process  A  consists  of  the  focal 
plane  coordinates  for  each  star  image.  Since  Process  A  calcu¬ 
lations  for  one  data  frame  can  be  performed  in  near  real  time 
and  many  times  faster  than  the  attitude  can  be  determined,  it 
will  likely  prove  desirable  to  perform  additional  editing  of 
the  star  data.  For  example,  images  with  rapidly  varying 
image  intensity  from  frame  to  frame  could  be  eliminated  or 
images  whose  successive  positions  are  inconsistent  with  the 
overall  motion  caused  by  vehicle  motion  (such  as  images  of 
space  debris)  could  be  deleted  immediately  from  consideration 
(failure  to  detect  and  delete  all  spurious  images  does 
not  prove  fatal,  but  does  slow  the  pattern  recognition 
logic  of  Process  B) .  Process  A  would  be  expected  to  output 
image  coordinate  and  magnitude  data  at  the  rate  of  about  5 
frames  per  sensor  each  minute  and  simply  overwrite  old  data. 
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The  microcomputer  of  Process  A  is  considered  as  an  integral 
part  of  the  star  tracker  itself,  making  it  a  "smart  sensor". 

Since  Process  A  controls  the  scan  of  the  CCD  and  its 
electronics,  it  is  possible  to  track  only  those  stars  desired 
(those  whose  pixel  response  lies  within  specified  bounds).  Thus, 
even  though  the  CCD  array  contains  thousands  of  pixels,  only  a 
small  fraction  of  their  response  values  need  be  subjected  to 
A/D  conversion  and  stored  at  any  one  time.  It  is  this  data 
compaction  feature,  along  with  the  fact  that  pixel  locations 
can  be  accurately  calibrated,  that  make  CCD  arrays  so  attractive 
for  this  application.  Also  significant  is  that  the  high  speed 
readout  of  the  CCD  allows  one  to  assume,  for  most  cases,  that 
the  star  images  visible  in  a  given  frame  have  been  imaged 
simultaneously.  It  is  therefore  possible  to  use  stellar 
resection  (geometric)  methods  and  not  necessary  to  account  for 
the  vehicle  motion  between  successive  "hits"  within  a  single 
frame  of  data. 

3)  Microcomputer  B 

Data  from  Process  A  (and  Process  C)  are  analyzed  by 
program  PtiOC&AA  8;  again  by  means  of  a  dedicated  microcom¬ 
puter.  As  input.  Process  B  accepts: 

*  star  image  coordinates  and  magnitude  data;  one 
set  per  FOV  (from  Process  A), 

*  a  priori  attitude  estimates  and  covariance,  (from 
Process  C) ,  and 

*  a  priori  estimates  and  covariance  of  interlock  angles 
between  the  sensors  image  planes  (from  Process  C) . 
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The  sequence  of  calculations/logical  decisions  divides 
into  two  primary  functions: 

*  identify  measured  stars  in  each  FOV  as  specific  stars 
contained  in  an  on-board  star  catalog  (containing,  in 
the  general  case,  the  direction  cosines  and  instrument 
magnitudes  of  the  5000  brightest  stars)  and 

*  determine  the  spacecraft  orientation  and  field  of  view 
interlock  angles  which  cause  the  simulated  images  of 
identified  catalog  stars  to  overlay  the  corresponding 
measured  images  in  a  least-squares  sense. 

These  two  tasks  will  be  discussed  in  detail  in  Sections 
2.1  and  5.1.  The  expected  output  rate  for  Process  B  is  two  or 
more  attitude  updates  per  minute  of  elapsed  time.  The  old 
attitude  and  covariance  is  overwritten  by  each  new  attitude 
and  covariance  and  the  new  values  are  output  to  Process  C. 

4)  Microcomputer  C 

The  attitude  determined  by  Process  B  for  a  discrete  time 
is  further  processed  by  program  PtiOC&AA  C  in  microcomputer  C. 
Input  to  this  program  consists  of  the  attitude  and  covariance 
from  Process  B  and  A/D  converted  gyro  rate  measurements  of 
angular  velocity.  The  previous  orientation  estimates  are 
integrated  forward  (using  the  gyro  rate  measurements)  to  the 
time  associated  with  the  new  orientation  from  Process  B.  The 
two  orientations  are  combined  in  a  Kalman  Filter  calculation 
to  give  a  best  estimate  of  real-time  attitude.  Further 
forward  integration  gives  an  estimated  orientation  and  covariance 
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for  the  next  Process  B  input. 

During  Phase  II  of  this  project  our  efforts  have  been  primarily  in 
the  following  areas: 

1)  Conversion  from  Euler  angles  to  Rodrigues  parameters 
to  describe  the  spacecraft  orientation  (affecting  the 
algorithms  of  Processes  B  and  C) . 

2)  Improvements  in  the  efficiency  of  Processes  B  and  C 
calculations  and  adding  several  needed  refinements. 

3)  Formulation  of  a  "truth  model"  for  adequately  testing 
Processes  B  and  C. 

4)  Design  of  an  efficient  star  catalog  format  for  easy 
access . 

5)  Develop  the  formulation  for  conversion  from  Rodrigues 
parameters  to  Euler  parameters  for  satellite  attitude 
in  both  Processes  B  and  C. 

These  and  other  topics  are  discussed  in  this  report.  All  our  work  involving 
computer  software  during  Phase  II  has  been  performed  using  the  HP9845A  micro¬ 
computer  system.  We  also  discuss  herein  the  anticipated  continuation  of  the 
above  work  and  new  tasks,  to  be  carried  out  during  Phase  III  of  the  project. 
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2.0  ORIENTATION  PARAMETERS 

The  Phase  I  report  outlined  the  Process  B  calculations  needed  to 
determine  the  vehicle  orientation  from  the  star  image  coordinates.  The 
solution  equations  involved  a  set  of  three  Euler  angles  defining  the 
rotation  sequence  for  orienting  the  field  of  view  (FOV)  relative  to 
inertial  space.  Any  3  angle  set  leads  to  differential  equations  and 
estimation  algorithms  containing  two  geometric  singularities  (<p  =  ±  n/2 
for  the  1-2-3  u-<p-K  sequence.)  Also,  Euler  angles  are  invariably  governed 
by  equations  having  trigonometric  nonlinearities;  this  increases  the 
expense  of  numerical  integration.  A  significant  change  in  Process  B 
since  the  Phase  I  report  has  been  the  introduction  of  Rodrigues  parameters 
as  a  set  of  3  parameters  defining  the  FOV  orientation.  The  Euler-Rodrigues 
parameters  have  only  one  singular  orientation  and  they  are  governed  by 
equations  having  quadratic  nonlinearities.  During  Phase  III  of  this 
project  we  will  explore  the  use  of  Euler  parameters  to  define  FOV  orien¬ 
tation.  The  four  Euler  parameters  have  the  advantages  that  (1)  they  do 
not  have  a  geometric  singularity,  and  (2)  they  rigorously  satisfy  linear 
differential  equations.  The  presence  of  the  constraint  that  they  sum 
square  to  unity  causes  some  problems  in  estimation  algorithms;  however, 
we  recently  found  a  method  to  circumvent  this  difficulty.  Since  Rodrigues 
parameters  are  most  easily  derived  from  Euler  parameters,  we  shall 
describe  the  latter  parameter  set  first. 

2.1  Euler  Parameters 

Euler  parameters  can  be  interpreted  geometrically 

in  terms  of  Euler's  theorem:  "A  completely  general  angular  displacement 
of  a  rigid  body  can  be  accomplished  by  a  single  rotation  (the  principal 
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angle,  <|>)  about  a  line  (the  principal  line,  l)  which  is  fixed  relative 

A  A  A 

to  both  arbitrary  body-fixed  axes  {b}  and  reference  axes  {n}.  If  {n} 
Is  initially  coincident  with  { b } ,  then  the  direction  cosines  (£j, 
of  n  with  respect  to  {n}  and  { b }  are  identical." 

The  Euler  parameters  are  then  related  to  the  principal  rotation 
parameters  as  follows: 


0Q  =  cos  <J>  / 2 

0^^  =  sin  <j>/2  ,  i=*l,  2, 3. 


(2.1) 


Note  that  Euler  parameter  satisfy  the  constraint: 
3  ? 

I  4  - 1. 

i=0 


(2.2) 


The  rotation  matrix  [C]  characterizing  the  relationship  between  a  body 
fixed  frame  { b }  and  a  reference  frame  (n)  by:  {b}  =  [C]{n}  can  be 
written  in  terms  of  Euler  parameters  as: 


[C]  = 


2  2  2  2 
B0  +  B1  ~  B2  "  B3 


2(B1B2  "  B0B3} 


2 (BiB3  +  W 


2(0102  +  0O03) 


2  2  2  2 
80  -  81  +  ®2  -  ®3 


2(82B3  "  “o8^ 


2(6363  -  e0e2) 


2(6363  +  6,363) 


2  2  2  2 
60  "  B1  "  B2  +  6 3 


(2.3) 


The  two  main  advantages  of  Euler  parameters  are  that  (1)  no  singularity 
exists  for  any  orientation,  and  (2)  no  evaluation  of  trignometric  functions 
need  be  done.  The  constraint  (2.2),  must  be  carefully  accounted  for, 
however . 

We  intend  to  explore  the  use  of  Euler  parameters  for  both  Process  B 
and  C  during  Phase  III.  From  preliminary  analysis  it  appears  Process  B 
can  be  adapted  easily  to  Euler  parameters.  The  only  disadvantage  in 
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using  Euler  parameters  is  the  need  for  a  constraint  equation  in  performing 
the  least-squares  correction.  We  outline  two  methods  we  have  tried  for 
incorporating  this  constraint  equation  in  Process  B. 

In  Process  B  we  seek  to  minimize  the  sum  of  the  squares  of  the 
residuals  between  measured  star  image  coordinates  and  predicted  coordinates 
for  the  same  stars,  using  direction  cosines  from  the  on-board  catalog. 

The  mapping  of  catalog  positions  onto  the  CCD  image  plane  is  a  function 
of  the  orientation  parameters  (Euler  angles,  Euler  parameters,  or  Rodrigues 
parameters)  via  the  stellar  colinearity  equations: 

milLl  +  MUL2  + 

A^31L1  +  AN32L2  +  AN33L3  ) 

(2.4) 

+  AV2  + 

y  '  [M31L1  *  "SzH  +  AH33L3  f 

where  f  a  lens  focal  length 

AN^j  =  elements  of  the  coordinate  frame  rotation  matrix  [AN]. 

*  star  direction  cosines. 

If  we  let: 

X  =  {(x^.y^)}  =  vector  of  calculated  CCD  image  plane  coordinates. 

X  *  {(x,,y,)}  =  vector  of  measured  star  image  positions  on  the  CCD. 

llm  - 

and  AX  =  X  -  X  ■  vector  of  residuals, 
then  we  seek  to  find  the  set  of  Euler  parameters,  B,  such  that  the 
weighted  sum  of  the  squares  of  the  residuals  is  minimized;  i.e.,  minimize 
<P  =  AXTWAX  (2.5) 

It  can  be  shown  that  minimizing  the  function  $  can  be  achieved  by  mini¬ 
mizing  a  function  $  iteratively;; 

<b  -  AXTWAX  (2.6) 

P  P  P 

where  AX  *  X  -  X  and 
P  P 
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Xp  *  vector  of  linearly  predicted  image  coordinates* 

But,  by  first-order  Taylor  expansion 

AX  -  AX  -  A  Ag 
P  c 

where 

AX^  **  vector  of  current  image  coordinates  residuals  based 
current  estimates  of  8* 

A  -  matrix  of  partial  derivatives  of  the  colinearity 
equations  with  respect  to  Euler  parameters 

A8  “  corrections  to  the  current  estimates  of  Euler  parameters. 

Thus,  we  can  write: 

$p  -  (AXc  -  AA8)T  W(AXc  -  AA8).  (2.7) 

In  addition  to  finding  the  set  of  Euler  parameters  to  minimize  ,  we 
must  also  satisfy  the  constraint  equation: 

0T0  -  1.  (2.8) 

Letting  8p  *  0c  +  AB,  we  find: 

(0  +A0)T  (0  +A8)  =*  1 

c  c 

or  to  first  order: 

1  '  BcBc  =  2BcAB-  (2-9) 

Thus,  our  problem  requires  that  we  minimize  Eq.  2.7  subject  to  the 

constraint  equation,  Eq.  2.9.  To  obtain  the  corrections  we  minimize  the 

augmented  function,  <t>,  using  the  method  of  Lagrange  multipliers: 

t  -  (AX  -  AAB)T  W(AX  -  AA0)  +  A (2AY  -  20TA0)  (2.10) 

c  c 

where  AY  -  (1  -  0T0)/2 

A  -  a  scalar  Lagrange  multiplier. 

For  minimization,  we  require: 

-  -2ATWAXc  +  2(ATWA)A0  -  2A0  -  0  (2.11) 

and 
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Vx$  -  2AY  -  26T  =  0  or 

T  (2.12) 

AY  -  8  AS. 

We  combine  these  two  equations  as  follows.  Solve  Eq.  2.11  for  A6: 

AB  -  (ATVA)-1ATWAXc  +  (ATWA)_1BX  (2.13) 

Substitute  this  into  equation  2.12: 

AY  -  BT(ATWA)-1ATWAXc  +  BT(ATWA)-1BA.  (2.14) 

Solve  for  A: 

-1  $ 

A  -  (6T(ATWA)_18)  (AY  -  6T(ATWA)_1ATWAX  ).  (2.15) 

c 

Finally,  substitute  for  A  in  Eq.  2.13  to  eliminate  A.  We  prefer  to 
write  this  final  equation  as: 

A0  =  A(f  +  K(AY  -  BTA8)  (2.16) 

where  A?  -  (ATWA)  ATWAX 

c 

=  corrections  to  B  without  the  constraint  and 
K  *  (ATWA)  B(BT(ATWA)  b)-1. 

Unfortunately,  Eq.  2.16  does  not  work  well  for  our  application,  although 

T 

the  general  formulation  is  useful  in  other  applications.  Matrix  A  WA 

appears  to  be  singular  or  very  nearly  so  and  thus,  in  the  present  appli- 
T 

cation,  (A  WA)  can  not  be  determined  (matrix  A  "knows"  that  there  are  really 
only  3  degrees  of  freedom!). 

Therefore,  we  have  tried  an  alternate,  and,  in  the  limit,  equivalent, 

solution  for  AB*  The  constraint  equation  can  be  incorporated  as  an  additional 

T 

"perfect"  observation  equation  but  with  a  large  weight.  That  is,  AY  *  (1  -0  0)/2 

T 

is  appended  to  the  AXc  vector  of  equation  (2.7)  and  B  appended  as  an  addi¬ 
tional  row  in  to  the  A  matrix  of  equation  (2.7).  The  weight  for  this  equation, 

.  T  "I 

the  (4,4)  element  of  W,  is  chosen  large  enough  (about  10® )  so  that  (A  WA)  does  not 
change  appreciably  for  variations  in  this  weight.  The  solution  for  AB  is 
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(2.17) 


then  the  classical  normal  equations: 

AB  -  (ATWA)  ATWAX 

c 

This  solution  for  AB  has  been  found  to  be  stable  for  all  cases  we  have 
tried  and  Is  more  convenient  than  equation  (2.16)  (which  is,  by  the  way, 
simply  th6  Kalman  Filter  algorithm  particularized  to  the  present  case  of 
a  perfect  measurement  having  infinite  weight  or  zero  variance). 

The  use  of  Euler  parameters  in  Process  C  should  be  a  straight¬ 
forward  modification  of  the  present  software,  which  involves  Rodrigues 
parameters.  In  addition,  there  are  several  properties  of  Euler  para¬ 
meters  which  will  be  utilized  to  advantage. 

We  have  decided,  for  the  present  formulation,  to  use  a  seven 
parameter  state  vector  for  Process  C.  This  consists  of  the  four  Euler 
parameters  (B^(t)}  and  three  unknown,  slowly  varying,  rate  gyro  biases, 
one  per  gyro.  The  biases  will  absorb  unmodeled  errors  such  as  gyro  non¬ 
orthogonality,  gyro  interlocks,  etc.  (Rigorous  justification  for  including 
only  the  rate  biases  is  difficult,  so  we  shall  do  simulations).  We  assume 
that  the  biases  are  constant  over  several  minutes  and  we  then  monitor  their 
values  over  longer  periods. 

The  Euler  parameter  differential  equations  describing  the  rotational 


kinematics  of  a  space  vehicle  are: 


1 

2 


Mt)] 


(2.18) 


f 


where  [w(t) ]  -  j 


-ui 


3 

u>2 


-w. 


Lw3 


where  {w^}  are  the  rate  gyro  readout  values  (including  noise)  for  vehicle 
rotation  with  components  along  the  3  orthogonal  gyro  frame  axes.  Since  the 
biases  are  assumed  constant  over  short  time  intervals  their  differential 
equations  are: 


This  set  of  seven  equations  is  integrated  forward  in  time  using  4-cycle 
Runge-Kutta  methods.  The  values  of  {{^(t)}  are  then  passed  to  Process  B 
to  be  used  for  the  analysis  of  the  next  frame. 

In  addition  to  this  set  of  equations,  we  must  also  obtain  the  (7x7) 
state  transition  matrix,  4>(t,t^),  for  the  Kalman  filter  state  update 
equations: 


[<Ht,tk)] 


3(6Q(t) ,  ei(t),  B2(t),  B3(t),  b1(t),  b2(t)  b3(t)) 
3(BQ(tk),  61(tk),  B2(tk),  S3(tk),  b2(tk),  b3(tk)) 


(2.20) 


The  upper  left  4x4  portion,  of  $  can  be  found  from  integration  of: 

t-  ~  dt— ]  -  [«o(t)J  [4»11(t, tk)J  .  (2.21) 
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(2.23) 


and  satisfies: 

(6<t))  -  l*n(t,tk)HfKtk)}  ,  *u(tk,tk)  »  [I). 

However,  detailed  analysis  of  the  16  differential  equations  comprising 
Eq.  2.2  reveals  that  there  are  only  4  distinct  equations  for  elements 
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(2.24) 


Equation  2.23  becomes,  then, 


this  can  be  rearranged  ("transmuted")  to: 


\ 

W 

w 

W 

B3(tk) 

l 

( 

W 

-w 

w 

-62Ctk> 

I 

II 

B2(tk> 

-w 

w 

J 

) 

-W 

e2(tk) 

-8o(tk>_ 

I 

(2.25) 


By  noticing  that  the  coefficient  matrix  above  is  orthogonal  we  can 
write  its  inverse  as  simply: 
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60(tk> 


ei(tk) 


62(tk) 


63(tk>' 


6l(tk> 


-e3(tk) 


B2(tk) 


02^tk^ 


W 


‘W 


-w 


63(tk) 


w 


B0(t) 


B:(t) 


B  2  ( t ) 


B3(t) 


Thus,  the  elements  of  *^2.^t,tk^  Can  °htained  from  the  values  of 


The  4x3  matrix  4>^2  must  be  formed  by  direct  integration  of 


i  t  ( t ,  t.  ) 

. — =  [u.(t)][«12(t,tk)]  + 


[Ai2] 


where 


fi(B0(t),  Bx(t),  e2(t),  B3(t)f| 
[A12]  =[3(bi,  b2,  b3)  J 


[$12(tk’tk)1  =  °* 


Kalman  Filter  Section  of  Process  C. 


The  state  covariance  and  Kalman  gain  update  recursions  are  the 
same  forms  used  in  the  Phase  I  report: 


Xk+1  “  X^k+l)  +  K(k+l){Y(k+l)  -  Y (k+1)  } 


Pk(k+1)  =  $(k+l,k)  Pk(k)<t  (k+1, k)  +  Q(k) 


K,.  ...  =  P,  (k+l)HT(k+l)[A  +  H(k+l)P  (k+l)HT(k+l)] 

U+i;  K  vk+lVk+l 


-1 


Pk+1(k+l)  =  [I  -  K(k+l)H(k+l)]Pk(k+l) 


H ( k+1 )  = 


3Y 


3X 


X^k+l) 


(2.26) 


{ 0 . ( tk> }  and  { B ^ ( t ) }  rather  than  from  integration  of  the  original  16 


equations  of  This  property  will  be  used  to  advantage  to  reduce 

the  number  of  calculations  in  Process  C. 


(2.27) 


(2.28) 


(2.29) 
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where  X^Ck+l) 


state  X  at  time  t^+^  based  upon  k  data  subsets 
forward  Runge-Kutta  integration  of  equations  2.18  and 
2.19  from  the  previous  estimates  X^(k), 


Y  (k+1)  = 


VW 

VW 

|92(tk+l> 

S3(tk+1) 


first  four  elements  of 


Xj^k+1). 


Y(k+1)  =  values  of  (8(t^+^)}  from  Process  B,  and 


Q(t)  =  process  noise  covariance  matrix  to  reflect  uncertainty  in 
(8q»  8^,  $2’  63)  arising  due  to  integration  of  errors  in 
the  rate  gyro  measurements. 


2.2  Rodrigues  Parameters 

In  constrast  to  Euler  parameters,  Rodrigues  parameters  form  a  three 
parameter  set,  thereby  eliminating  the  need  for  a  constraint  equation. 
However,  the  Rodrigues  parameters  have  a  singularity  like  Euler  angles 
(although  any  Euler  angle  set  has  two!).  Computationally  they  are 
attractive  since  they,  too,  do  not  involve  evaluation  of  trignometric 
functions  for  either  the  least-squares  correction  of  Process  B  nor  the 
Runge-Kutta  integration  of  Process  C. 

Rodrigues  parameters,  {q^}  are  related  to  Euler  parameters  {8^} 
by: 

qi  =  V6o  =  *i  tan  41/2  ’  i=1*2’3  (2.30) 


The  singularity  occurs  at  <t>  =  180°.  The  3*3  matrix  for  coordinate  frame 
rotation,  [C],  expressed  in  Rodrigues  parameters  is: 


[C] 


2(qxq2  -  q3) 
2(qjq3  +  q2) 


2(qxq2  +  q3) 

i  2  .  2  2 

1  "  ql  q2  ~  q3 
2 (q2q3  "  ql) 


2(q1q3  -  q2) 

2(q2q3  +  q3) 

.  2  2  2 
1  _  ql  _  q2  +  q3 


(2.31) 


and  the  kinematic  equations  are: 
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lw3  / 

(2.32) 


Process  B  and  C  are  currently  operational  using  these  parameters. 

The  simulations  described  in  this  report  have  all  been  performed 
using  the  Rodrigues  parameters.  Process  B  and  C  equations  are  very  similar  to 
the  development  reported  in  section  2.1  except  for  3  parameters  instead 
of  4.  The  tradeoffs  between  Euler  parameters  and  Rodrigues  parameters 
are  the  following: 

(1)  The  Euler  parameters  are  universal;  no  geometric 
singularity  exists ,  whereas  the  Rodrigues  parameters 
still  possess  a  singular  point. 

(2)  The  Euler  parameters'  state  transition  matrix  can  be 
found  without  numerical  integration  [see  Eqs.  (2.21)  - 
(2.26)],  whereas  these  developments  do  not  hold  for 
Rodrigues  parameters. 

(3)  The  three  Rodrigues  parameters  result  in  lower  dimen¬ 
sional  matrices  to  invert  than  the  four  Euler  parameters, 
in  Process  B  calculations. 


Advantages  (1)  and  (2)  outweigh,  in  our  judgement,  disadvantage  (3)  of  the 
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1 


Euler  parameters  in  comparison  to  Rodrigues  parameters.  It  is  also 
significant  to  note  that  both  Euler  parameters  or  Rodrigues  parameters 
are  computationally  superior  to  Euler ian  angles,  since  no  trigonometric 
nonlinearities  are  present  in  the  differential  equations. 
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3.0  COORDINATE  FRAMES 


We  have  used  several  coordinate  frames  thus  far  for  the  development 
of  Process  B  and  C  algorithms.  Refer  to  the  Phase  I  Report  for  a  complete 
discussion  of  possible  frames  and  their  interrelationships. 

The  two  FOV  frames  used  for  the  research  described  here  are  assumed 
to  be  separated  by  90°  and  their  orientations  are  nominally  fixed  relative 
to  the  vehicle  gyro  axes.  The  gyro  frame  unit  vectors,  {g^},  are  assumed 
oriented  such  that  is  along  the  vehicle  x-axis  (pointing  in  the  direction 
of  flight)  and  g2  alK*  £.3  are  perpendicular  to  the  vehicle  x-axis.  The 
FOV(A)  frame  orientation  can  be  described  by  a  (1,2,3)  Euler  angle 
sequence  of  (90°,  135°,  0°)  relative  to  the  gyro  frame.  This  sequence 
places  a_2  (the  image  plane  "y"  axis)  along  ant*  £3  (the  boresight) 
rotated  45°  from  g^  about  g^.  The  FOV(B)  frame  is  oriented  relative  to 
FOV (A)  by  a  (3,1,3)  Euler  angle  sequence  of  (90°,  -90°,  90°).  Again, 
b^  lies  along  £  and  b^  (the  boresight)  is  rotated  45°  from  £  about  g^ 

(see  Figure  3.1). 

The  gyro  frame  to  FOV (A)  and  (B)  rotations  described  above  have 
been  used  for  generating  simulated  data  for  Process  B  (see  Chapter  4). 

The  gyro  frame  is  oriented  such  that  g2  *s  perpendicular  to  the  orbit 
plane  and  g^  along  the  radius  vector.  In  Process  B  we  do  not  use  the 
gyro  frame.  Rather,  the  algorithms  have  been  written  to  solve  for  the 
orientation  parameters  relating  FOV (A)  and  (B)  to  the  inertial  frame. 

We  have  not  pursued  the  use  of  an  orbiting  frame  as  suggested  in 
the  Phase  I  report.  The  expected  conversion  to  Euler  parameters  will 
obviate  the  need  for  defining  the  vehicle  orientation  relative  to  an 
orbiting  frame  because  Euler  parameters  contain  no  inherent  singularity 
as  do  Euler  angles  and  Rodrigues  parameters.  Thus,  by  using  Euler  para- 
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choice  of  a  suitable  orbiting  frame  to  avoid  singularities  for  each  type 


S3  ’-2  ’-2 


Figure  3.1.  Relationship  of  the  gyro  frame  to  FOV(A)  and  FOV(B) 


4.0  TRUTH  MODELS 


4.L  Process  B  Models 

In  order  to  adequately  test  Process  B  we  must  generate  realistic 
input  data  (i.e.,  realistic  output  data  from  Process  A).  This  requires 
generating  a  sequence  of  data  frames  using  an  assumed  vehicle  orbit  and 
appropriate  time  intervals  between  successive  frames.  The  image  coordinates 
must  be  perturbed  by  relevant  noise  sources,  and  occasional  additions  of 
spurious  images  as  well. 

The  first  step  in  the  data  generation  is  to  choose  an  appropriate 
spacecraft  orbit  by  selecting  the  semi-major  axis  of  the  orbit,  orbital 
period  and  initial  position  in  inertial  space.  To  facilitate  the  cal¬ 
culation  of  subsequent  satellite  positions  we  make  use  of  Herrick's  "f 
and  g"  solution  (see  Junkins,  p.  155).  This  method  is  used  also  for 
computing  the  position  and  velocity  of  the  earth  at  subsequent  times. 

Tests  to  date  have  been  for  a  nearly  circular  orbit  using  a 
nominally  earth  pointing  spacecraft.  The  nominal  gyro  frame  orientation 
is  then  determined  from  the  vehicle  velocity  and  coordinate  vectors. 

Thus,  is  along  \r,  ^  is  perpendicular  to  the  orbit  plane  and  will 
be  defined  by  the  cross  product  of  and  The  FOV (A)  and  (B)  frames 

are  then  defined  relative  to  {g}  as  described  in  section  3. 

Given  the  boresight  of  each  FOV,  the  star  catalog  is  accessed  in 
the  same  manner  as  in  Process  B  (see  Appendix  A  and  the  Phase  I  report). 

The  combined  velocity  of  the  earth  and  vehicle  is  used  to  compute  the 
aberration  of  starlight  for  each  star  (see  Appendix  B) .  The  resulting  star 
direction  cosines  are  then  rotated  into  the  FOV  and  the  image  plane  coordinates 
are  computed  using  the  colinearity  equations  and  an  assumed  lens  focal 


length  of  70mm  +  e,  where  e  is  a  small  defocuslng  length  needed  co  blur 
the  star- image  enough  to  allow  accurate  centroid ing.  For  our  purpose, 
we  have  assumed  that  the  centroiding  can  be  performed  to  an  accuracy  of 
better  than  10%  of  a  pixel  and  corrections  to  image  location  can  be 
applied  for  causes  such  as  image  distortions.  Thus,  we  add  random 
Gaussian  noise  (with  o  =  3.4u)  to  the  image  coordinates  of  each  star  to 
simulate  centroiding  errors.  Star  magnitudes  are  obtained  from  the  star 
catalogue  and  are  perturbed  with  noise  also. 

The  true  orientation  and  image  coordinates  as  well  as  the  error- 
corrupted  simulated  data  are  stored  on  computer  tape.  At  present,  we 
separate  data  frames  by  30  seconds  of  satellite  motion  since  that  is  the 
time  needed  for  Process  B  solution.  Upon  application  of  the  pattern  recog¬ 
nition  and  estimation  algorithms,  all  orientation  parameters,  interlock 
parameters,  etc.,  can  be  compared  versus  the  stored  truth  model. 

4.2  Process  C  Models 

For  our  tests  of  Process  C  we  generated  data  for  a  series  of  satellite 
positions,  again  separated  by  30  seconds  of  real  time.  The  data  con¬ 
sists  of  true  orientation  parameters  and  perturbed  values  representing 
output  from  Process  B  and  an  associated  covariance  matrix  which  would 

result  from  the  least-squares  solution.  Eventually,  we  will  use  the 

* 

actual  output  from  Process  B  as  input  to  Process  C.  Rather  than 
generating  and  storing  a  large  number  of  gyro  rate  samples  on  tape,  we 
instead  compute  the  rates  at  each  step  of  the  Runge-Kutta  integration 
and  add  Gaussian  noise  to  the  true  values,  to  simulate  gyro  and  A/D 
conversion  errors. 

Although  our  "truth  model"  for  Process  C  is  not  sophisticated,  we 

*The  linkup  of  Processes  B  and  C  can  be  implemented  within  the  same  micro¬ 
computer  or  using  parallel  microcomputers.  We  will  address  this  issue  in 
the  Phase  III  effort. 
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nevertheless  believe  it  is  satisfactory  for  testing  the  validity  of 
Process  C.  The  reason  for  this  belief  is  that,  in  a  steady  state  mode, 
the  integrated  state  will  usually  be  much  more  poorly  determined  than 
the  "fresh"  data  from  Process  B.  Thus,  the  optimal,  up-dated  state 
output  from  Process  C  will  nearly  equal  the  Process  B  output  data.  This 
is  borne  out  by  our  simulations. 
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5.0  SIMULATIONS 


5.1  Process  B 

The  functions  of  Process  B  are  outlined  in  Figure  5.1.  Presently, 
input  data  for  FOV(A)  is  read  from  tape;  this  consists  of  estimated 
orientation  parameters,  number  of  stars  in  each  FOV  and  a  list  of  image 
plane  coordinates  and  intensities  (magnitudes)  for  the  measured  stars. 

Given  the  orientation  estimates,  the  star  catalog  is  accessed  to  retrieve 
those  stars  which  lie  within  some  uncertainty  border  around  the  estimated 
boresight.  Cosines  of  inter-star  angles  are  computed  for  the  measured 
stars  and  tabulated.  We  next  begin  a  series  of  tests  comparing  cosines 
of  inter-star  angles  for  pairs  of  catalog  stars  with  the  table  of  measured 
values  and  look  for  a  match.  This  pairing  begins  with  stars  nearest  the 
estimated  boresight  and  proceeds  outward  until  a  match  is  obtained. 

Refer  to  the  Phase  I  report  for  more  details. 

Once  a  match  has  been  obtained  we  perform  a  least-squares  differential 
correction  to  correct  the  assumed  orientation  parameters  until  the 
mathematical  projection  of  identified  catalog  stars  onto  the  image  plane 
overlays  the  measured  stars  in  a  least-squares  sense.  At  this  point,  the 
subset  of  stars  from  the  catalog  is  searched  for  additional  confirming 
stars  (i.e.,  stars  that  project  onto  the  image  plane  near  to  a  measured  star). 
Stellar  aberration  effects  are  added  to  the  catalog  stars  and  a  final  least- 
squares  correction  is  performed  with  up  to  5  stars.  The  above  process 
is  repeated  for  FOV(B)  in  like  manner.  If  we  assume  the  interlock 
angles  between  FOV (A)  and  (B)  are  fixed  or  slowly  varying,  it  is  possible 
to  estimate  the  orientation  to  within  5  arc  seconds  for  the  errors 
assumed  in  the  image  coordinates  (o  '  3.4u). 


Figure  5.1:  OUTLINE  OF  PROCESS  B  FUNCTIONS 


FOV(Aj 

.  -  Read  data  from  buffers: 

-  measured  coordinates  from  Process  A 

-  orientation  estimates  from  Process  C 

-  Compute  (3x3)  rotation  matrix 

-  Obtain  estimated  boresight  vector  from  last  row  of  rotation  matrix 

-  Access  the  star  catalog  for  stars  surrounding  the  boresight 

(a  sub-catalog) 

-  For  measured  stars: 

-  calculate  cosine  of  interstar  angle  between  each  star  pair 

-  store  cosine  values  and  star  pair  indicator  in  a  table 

-  For  sub-catalog  stars: 

-  calculate  cosine  between  star  pairs 

-  compare  each  cosine  with  the  table  of  measured  cosines 

-  test  for  a  match  (if  no  match,  restart) 

-  Perform  iterative  least-squares  differential  correcton: 

-  calculate  the  matrix  of  derivatives,  (A),  consisting  of  partials 

of  stellar  colinearity  equations  with  respect  to  the  orientation 
parameters,  evaluated  for  the  two  paired  stars 

-  calculate  residuals  (AY)  between  measured  coordinates  and 

coordinates  for  the  two  catalog  stars  projected  on  the 
image  plane 

T  T 

-  compute  corrections  to  orientation  parameters:  Af?  =  (A  WA)  A  WAY 

where  W  is  a  weight  matrix 

-  add  corrections  to  current  values  =  6^  +  A6 

-  test  for  small  corrections;  repeat  if  necessary 

-  Search  for  confirming  star  images: 

-  compute  image  coordinates  for  all  stars  in  sub-catalog 

-  compare  coordinates  of  each  catalog  star  with  coordinates 

of  each  measured  star 


-  save  measured  and  catalogued  pair  if  close 

-  if  there  are  no  new  pairs  beyond  the  original  pair,  reject 

orientation  and  continue  cosine  pairing  to  look  for  another 
pair  match 

-  repeat  least-squares  with  up  to  5  stars 


FOV(B) 

-  Read  data  from  buffer: 

-  measured  coordinates  from  Process  A 

-  use  orientation  from  FOV(A)  solution 

-  Repeat  same  steps  as  for  FOV(A)  for  pairing  stars  and  least-squares 

correction 

Combined  Solution: 

-  compute  or  obtain  interlock  angles  between  FOV(A)  and  FOV(B) 

-  compute  derivative  matrix  (A)  for  up  to  5  stars  from  each  FOV 

-  compute  coordinate  deviations  (AY)  for  these  stars 

-  compute  correction  to  orientation  parameters: 

A6  =  (ATWA)  ATWAY 

-  add  corrections  and  test  for  smallness;  repeat  if  necessary 
Output: 

-  Output  to  buffer: 

-  orientation  parameters 

-  covariance  matrix  for  orientation  parameters 

-  associated  time 
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One  goal  of  the  Process  B  simulations  was  to  measure  the  average 
solution  time  for  a  set  of  data  from  Process  A.  Our  present  software, 
written  in  high-speed  interpreter  BASIC,  can  process  a  set  of  measurements 
in  about  30  seconds  of  real  time.  This  figure  includes  tape  and  disk 
read  time  (approximately  5  sec.).  The  total  of  30  seconds  is  certainly 
an  upper  limit  for  a  flight  system. 

A  more  important  goal  of  Process  B  simulation  is  to  measure  the 
orientation  accuracy  relative  to  input  errors.  Simulations  discussed 
here  involved  Rodrigues  parameters  which  do  not  yield  a  readily  interpre¬ 
table  result  for  error  analysis.  Therefore,  two  angles  were  computed 
which  give  a  clearer  measure  of  errors.  Since  the  true  orientation  is 
known,  we  can  compute  the  angular  error  in  boresight  direction  by  computing 
the  cross  product  between  the  true  and  the  calculated  boresights. 

This  angle  is  then  used  to  rotate  the  calculated  frame  into  the 
true  frame,  so  that  the  boresights  coincide.  Then  the  rotation  about  the 
true  boresight  is  calculated,  again  by  using  a  cross  product.  Refer  to 
Appendix  C  for  a  complete  discussion  of  this  procedure.  The  two  angles 
measure  the  angular  error  off  the  true  boresight  and  the  angular  rotation 
error  about  the  true  boresight.  The  results  of  several  typical  successive 
trials  are  listed  in  Table  5.1.  It  appears  from  this  table  that  5-10  arc 
second  accuracy  is  obtainable  for  the  assumed  geometry  and  system  parameters. 
As  will  be  discussed  more  fully  in  Section  6,  the  assumed  error  of  3.4w 
(or  1/10  of  a  pixel)  for  the  image  coordinates  is  probably  an  upper  limit 
and  may  be  reduced  to  5%  pixel  or  less  by  various  means.  Also,  the  error  is 
proportional  to  scale  factor,  i.e.,  attitude  error  is  inversely  proportional 
to  the  lens  focal  length  (and  field  size). 
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5.2  Process  C  Simulations 

The  software  for  Process  C  is  more  straightforward  than  that  of 
Process  B.  Process  C  is  functioning  presently  with  Rodrigues  parameters. 

The  state  vector,  then,  contains  six  elements:  3  orientation  parameters 
and  3  rate  gyro  bias  values.  A  typical  simulation  result  is  shown  in 
Figure  5-2.  To  test  Process  C  we  have  generated  simulated  output  from 
Process  B  consisting  of  three  Rodrigues  parameters  and  associated  covariance 
matrix  for  each  time  step  of  30  seconds.  Any  initial  set  of  three  Rodrigues 
parameters  is  suitable  for  testing  Process  C.  However,  generating  a 
realistic  covariance  matrix  associated  with  that  set  and  subsequent  sets  is 
difficult.  This  is  due,  in  part,  to  the  infinite  range  of  Rodrigues  para¬ 
meters  which  causes  a  wide  range  of  values  for  the  covariance  matrix.  Also, 
errors  in  orienting  the  FOV  are  not  easily  mapped  into  errors  in  the  Rodrigues 
parameters  and  the  covariance  matrix.  We  have  not  generated  rigorous 
covariance  matrices,  since  the  present  software  is  considered  a  "stepping- 
stone"  toward  the  Euler  parameter-based  software.  Therefore,  Figure  5-2 
should  be  viewed  as  illustrating  the  general  behavior  of  the  Kalman  filter 
in  response  to  data  from  Process  B  and  to  the  state  obtained  from  integrating 
the  rate  gyro  data. 

The  vertical  scale  of  Figure  5-2  is  in  arc-seconds,  approximately 
transformed  from  the  Rodrigues  parameters  output  from  the  Kalman  filter. 

We  have  also  included  (Figure  5-3)  Phase  I  results  of  Kalman  filter  simula¬ 
tions  using  Euler  angle  data.  From  a  comparison  of  Figure  5-2  and  5-3  one 
can  easily  see  the  similar  behavior  for  the  Kalman  filter  processing  of 
either  Euler  angle  or  Rodrigues  parameter  data. 

Our  simulations  have  shown  that,  after  a  steady  state  has  been  reached, 
the  updated  orientation  parameters  output  from  the  Kalman  filter  are  nearly 
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equal  to  the  parameters  output  by  Process  B.  The  near  equality  is  due 
to  the  significantly  smaller  covariance  associated  with  the  Process  B  para¬ 
meters  as  compared  to  the  covariance  associated  with  the  orientation  parameters 
obtained  by  forward  integration  of  the  kinematic  differential  equations  using 
rate  gyro  data.  However,  we  find  that  the  gyro  bias  estimates  obtained  from 
Process  C  oscillate  about  their  true  value  on  a  short  time  scale  and  with  a 
deviation  similar  to  the  gyro  noise  value  (refer  to  Table  5.2  and  Figure  5.2). 
This  oscillation  appears  to  be  due  primarily  to  the  similar  magnitudes  of  the 
gyro  biases  and  assumed  gyro  noise.  Thus,  gyro  noise  is  "interpreted"  bv  the 
Kalman  filter  as  a  change  in  the  gyro  biases  even  though  the  true  biases,  in 
these  simulations,  are  constant.  We  suggest  smoothing  the  raw  gyro  rates 
(by  averaging  multiple  values,  for  example),  since  the  gyro  readout  rate 
is  much  higher  than  required  by  the  Runge-Kutta  integration.  Alternatively, 
the  updated  bias  values  can  be  ,'veraged  into  a  running  average  (an 
average  of,  say,  the  last  50  values).  Application  of  either  or  both  of 
these  methods  would  smooth  the  recovered  bias  values  sufficiently  so  that 
we  can  monitor  long  term  bias  changes.  These  methods  have  been  applied 
for  short  series  of  runs.  During  Phase  III  we  will  perform  more  extensive 
tests. 
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Figure  5.3b  Kalman  filter  results  for  Process  C  simulation:  Case  II 
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Table  5.2  Recovered  bias  values  for  10  min  series. 
Initial  estimates  were  0. 


time (sec ) 

bl 

b2 

b3 

30 

5.761E-6 

8.822E-6 

6 . 694E-6 

60 

5.220 

7.076 

4.816 

90 

6.635 

2.755 

3.025 

120 

6.540 

3.802 

4.965 

150 

4.769 

-0.739 

5.057 

180 

5.412 

2.053 

4.646 

210 

5.449 

2.991 

4.180 

240 

5.  761 

4.873 

3.723 

270 

5.935 

5.200 

3.634 

300 

6.000 

5.311 

3.597 

330 

5.626 

5.285 

3.734 

360 

5.335 

4.935 

4.406 

390 

5.098 

4.424 

4.752 

420 

4.691 

3.959 

4.904 

450 

4.742 

3.775 

4.881 

480 

4.751 

3.752 

5.012 

510 

5.008 

3.798 

5.357 

540 

4.921 

3.913 

5.127 

570 

5.018 

4.086 

5.442 

600 

5.109 

4.219 

5.408 

Average 

5.389 

4.215 

5.151 

True  Value 


4.848E-6 


4.848E-6 


4.848E-6 


6.0  CONCLUSIONS  OF  PHASE  II  RESEARCH 

We  have  met  several  research  objectives  through  the  efforts  during 
Phase  II.  The  research  has  also  suggested  new  directions  for  future 
refinements  and  improvements. 

Our  simulations  of  Process  B  have  focused  on  "steady-state"  operation 
That  is,  simulations  in  which  the  estimated  orientations  were  not  in 
error  by  a  large  amount  and  there  were  at  least  several  stars  in  each 
FOV.  In  nearly  all  steady-state  examples  the  solution  time  was  30 
seconds  or  less  of  real  time,  including  about  5  seconds  required  for 
reading  input  data  from  tape  and  for  accessing  the  star  catalogue  stored 
on  disk.  The  solution  time  for  Process  B,  as  it  is  currently  configured, 
is  an  upper  limit  and  can  be  shortened  considerably.  There  are  several 
factors  leading  to  a  decrease. 

We  note  that  memory-to-memory  data  transfers  are  significantly 
faster  than  memory-to-storage  device  transfers.  Therefore,  a  flight 
computer  system,  employing  memory  buffers  for  data  transfers  between 
subsystems  and  a  memory  bank  for  the  star  catalog,  would  save  at  least 
five  seconds  over  the  present  system  over  and  above  savings  which  may  be 
realized  as  a  result  of  the  particular  flight  computer  doing  arithmetic 
faster  than  the  HP9845A. 

A  more  significant  time  savings  can  be  realized  by  programming 
Process  B  in  microcode  or  by  using  a  compiled  program  form.  The  present 
computer  system  uses  a  high-speed  BASIC  interpreter;  i.e.,  each  line  of 
BASIC  code  must  be  translated  or  interpreted  and  then  executed.  A 
factor  of  five  increase  in  speed  may  be  possible  with  a  microcode  program. 
On  the  basis  of  our  simulations  and  the  above  comments,  we  conclude  that 
Process  B  can  easily  be  carried  out  within  the  constraints  of  on-board 
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computations .  We  note  that  Process  B  (as  refined  during  Phase  III)  will 
have  some  additional  logic  and  this  will  lengthen  the  program  slightly. 

The  bulk  of  the  computations  are  being  performed  in  the  present  version, 
however. 

We  are  also  near  our  accuracy  goal  of  5-10  arcseconds  using  an  adopted 
system  and  associated  parameters  when  there  are  4-5  stars  per  FOV.  Before 
we  can  unequivocally  make  this  claim  for  accuracy,  we  must  be  sure  we  are 
accurately  modeling  Process  A  output  data  (star  image  centroids).  We 
emphasize  that  the  accuracy  of  the  Process  B  output  data  depends  directly 
on  the  quality  of  Process  A  output  data.  Thus,  any  accuracy  claim  for 
Process  B  must  be  accompanied  by  figures  for  the  accuracy  of  Process  A 
(and  type  of  star  tracker).  The  5-10  arcseconds  goal  can,  no  doubt,  be 
exceeded  by  redesigning  the  system  around  a  longer  focal  length  lens  and 
admitting  fainter  stars.  However,  the  storage  and  manipulation  of  the 
resulting  increased  population  of  visible  stars,  may  be  beyond  the  capa¬ 
bilities  of  a  real-time  on-board  system. 

Process  C  has  not  been  tested  as  extensively  as  has  Process  B.  How¬ 
ever,  our  simulations  of  Process  C  show  that  it  satisfies  the  time  con¬ 
straints  of  on-board  computation.  Process  C  is  extremely  fast  —  typically 
5  seconds  for  the  Runge-Kutta  integration  and  Kalman  filter  update  sections 
combined.  The  decrease  in  solution  time  by  using  a  microcode  version  will 
not  be  as  dramatic  as  for  Process  B.  Much  of  the  current  software  makes 
extensive  use  of  the  BASIC  matrix  commands  of  the  HP9845A  (such  as  matrix 
multiplication  and  matrix  inversion),  each  of  which  is  already  in  machine 
language  and  therefore  very  fast.  However,  we  expect  some  savings  by 
converting  to  Euler  parameters  and  using  the  associated,  more  efficient, 
relationships  for  computing  the  state  transition  matrix. 
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7.0  OUTLOOK  FOR  PHASE  III  RESEARCH 


Phase  III  research  will  focus  primarily  on  extensive  testing  and 
integration  of  the  complete  attitude  system.  We  will  continue  to  study 
new  approaches  and  new  techniques  as  well. 

Our  first  task  will  be  the  conversion  from  Rodrigues  parameters  to 
Euler  parameters,  as  outlined  in  Section  2.1.  This  task  should  be 
straightforward  since  most  of  the  calculations  and  logic  have  been 
formulated.  Software  conversion  has  begun. 

A  second  significant  task  is  the  editing  and  transfer  of  the  SKYMAP 
star  catalog  to  the  HP9845  system.  The  catalog  will  be  stored  on 
disk  in  place  of  the  present  simulated  catalog.  The  SKYMAP  catalog, 
currently  stored  on  the  IBM  370/158  computer,  may  be  modified  prior  to 
transfer.  No  attempt  was  made  originally  to  have  the  catalog  uniform 
in  number  of  stars  per  cell.  We  simply  selected  the  5000  brightest 
stars.  However,  for  accurate  attitude  determination  we  may  desire  at 
least  5  detectable  stars  per  FOV.  Near  the  galactic  pole  we  will  need  a 
fainter  limiting  magnitude  in  order  to  obtain  this  number  compared  to 
the  galactic  equator  regions.  To  maximize  the  overlap  between  the 
observed  stars  and  catalogued  stars,  the  integration  time  of  the  CCD 
sensor  must  be  variable  and/or  Process  A  logic  structured  to  select  and 
process  images  from  the  5  or  6  brightest  stars  in  any  field. 

We  must  also  delete  double  stars  from  the  catalog.  This  includes 
all  star  pairs  with  separations  less  than  the  resolution  of  the  CCD 
sensor.  We  have  postponed  this  task  until  we  have  a  firmer  resolution 
value  for  a  typical  CCD  star  sensor. 

Further  refinements  will  be  added  to  Process  B.  We  must  add  calcu¬ 
lations  and  logic  to  monitor  the  FOV  interlock  angles.  A  third  FOV 
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should  be  added,  at  least  for  general  use,  so  that  any  pair  of  sensors 
could  be  used  for  attitude  determination.  A  third  sensor  would  provide 
a  backup  should  one  sensor  fail  and  permit  a  small  sun  angle  for  one 
sensor,  leaving  two  others  operational  for  accurate  attitude  determination. 

Appropriate  "start-up"  logic  must  be  incorporated  into  Process  B; 
this  logic  decides  whether  to  continue  searching  for  a  star  pair  match 
or  to  reject  Process  A  data  and  request  a  new  set.  We  are  presently 
studying  techniques  with  which  to  re-access  the  star  catalog  should 
Process  B  fail  to  find  a  star  pair  match  near  the  estimated  boresight. 

Our  current  thoughts  are  to  access,  in  succession,  five  or  six  adjacent 
FOV-size  star  subsets  located  around  the  original  estimated  FOV.  Each 
star  subset  would  be  treated  in  the  same  manner  as  the  subset  containing 
the  original  estimated  FOV.  The  software  for  this  scheme  has  been  written 
but  not  yet  incorporated  into  the  existing  Process  B  code. 

As  our  final  task,  we  propose  to  run  extensive  simulation  tests  of 
Processes  B  and  C.  To  do  this,  we  will  first  simulate  as  accurately  as 
possible  the  data  from  Process  A.  As  mentioned  before,  some  of  the 
error  sources  can  be  corrected  for  in  Process  A.  There  will  be  errors 
in  positions  due  to  random  noise  inherent  in  the  CCD  and  associated 
electronics  and  noise  due  to  response  variations  from  pixel  to  pixel. 

Errors  of  this  nature  depend  on  the  type  of  CCD  and  electronic  design 
used.  They  are  presently  modeled  as  errors  in  star  position  which  belong 
to  Gaussian  distributions.  The  accuracy  results  of  Process  B  and  C 
-Hron1 »tions  will  then  be  proportional  to  the  selected  standard  deviation. 

We  will  assume  that  systematic  errors  in  the  star  sensors  can  be  removed 
by  apriori  and  dynamic  calibrations  applied  in  Process  A. 

The  simulated  data  from  Process  A  will  represent  data  frames  separated 
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by  perhaps  30  seconds  of  time  and  for  several  fractional  orbits  of  a 
typical  satellite.  Tests  will  determine  1)  the  rate  of  failure  for 
Process  B  to  determine  any  attitude  and  2)  the  accuracy  of  the  solutions 
for  attitude  as  a  function  of  number  of  stars,  possibly  the  placement  of 
the  stars  in  the  field,  and  magnitudes  of  the  stars.  By  randomly  pointing 
the  boresight  axes  we  can  also  test  the  "start-up"  capability  of  Process 
B  relative  to  initial  estimates. 

Process  C  tests  will  determine  the  ability  of  the  Kalman  filter  to 
produce  an  optimal  estimate  near  the  true  state  vector  in  the  steady 
state  mode.  Finally,  the  speed  of  Process  B  and  C  will  be  measured. 

We  shall  also  determine  whether  Processes  B  and  C  actually  need  to 
be  run  in  parallel  or  if  they  can  function  adequately  in  a  sequential 
fashion.  If  a  parallel  mode  is  desirable,  we  will  configure  a  parallel 
test /demonstration  of  Processes  B  and  C. 
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APPENDIX  A:  STAR  CATALOG  FORMAT 


As  of  this  writing  we  have  not  transferee!  the  SKYMAP  star  catalog 
from  the  IBM  370/158  computer  to  the  HP9845.  To  procede  with  testing 
and  to  gain  experience  with  the  form  for  the  catalog,  a  pseudo-catalog 
was  created  on  the  HP9845.  This  was  generated  by  computing  5000  random 
star  positions,  from  a  uniform  distribution  over  the  celestial  sphere, 
and  computing  star  magnitudes.  The  stars  were  then  reordered  into  cells 
using  the  techniques  described  in  the  Phase  I  report.  After  some  thought, 
it  was  decided  to  have  529  cells  (22  latitude  bands  of  equal  width) . 

Choice  of  this  number  was  influenced,  in  part,  by  the  FOV  size  of  7°  x 
9°  adopted  for  these  simulations.  With  this  field  size,  the  catalog 
access  routine  reads  data  from  the  4  nearest  neighboring  cells  around 
the  estimated  boresight  (Figure  A.l)  and  thus  provides  nearly  complete 
coverage  of  the  estimated  FOV  by  the  4  cells. 

To  facilitate  computer  access,  the  catalog  was  stored  on  disk  in 
a  direct  access  format.  Preceeding  the  catalog  proper  is  a  table  con¬ 
taining  the  position  in  the  file  of  each  cell  and  the  number  of  stars  in 
each.  Thus,  each  cell  can  be  accessed  rapidly  and  independently. 
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APPENDIX  B:  STELLAR  ABERRATION 


The  effect  of  stellar  aberration  is  to  cause  a  star's  apparent 
direction  to  shift  cowards  the  direction  of  the  observers  motion.  The 
amount  of  shift  depends  on  both  the  velocity  of  the  observer  and  on  the 
angle  between  the  observer's  line  of  sight  (the  star  direction)  and  the 
velocity  vector.  The  shift  is: 

a  *  -  sin  a  (B.l) 

c 

where  a  *  aberration  in  radians 
v  =  observer's  speed 
c  =  velocity  of  light 

a  =  angle  between  velocity  vector  and  the  true  star  direction. 


For  our  purpose,  we  must  express  a  star's  shifted  direction  in  terms  of 
the  true  direction,  the  vehicle  velocity  and  the  angle  between  the 
velocity  vector  and  true  direction;  that  is. 


x 

i 

•J 

y 

» 


=  f (Lx.Ly.Lz,v,a) 


If  we  let  v^  be  the  velocity  of  starlight  in  the 
v  be  observer  velocity,  then  the  relative  velocity  of 
seen  by  the  observer  is: 


(B.  2) 


inertial  frame  and 
the  starlight  as 


v  ,  =  v  -  v 

— s/o  —s  — 


(B.3) 


Now,  if  we  let  i  be  the  unit  vector  in  the  true  direction  and  i  the 
— n  — n 

unit  vector  in  the  shifted  direction,  we  can  rewrite  this  as  (refer  to  Figure 

Bl) 


(c  +  v  cos  a)  £'  =  cfc  +  v 
— n  -n  — 


(B.4) 
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To  first  order  this  becomes 

-  „ 

i '  =  (1  -  —  cos  ot)!l  +  — —  (K.5) 

— n  c  — n  c 

or 


This  equation  is  used  for  calculating  the  displacement  of  a  star's 
unit  vector.  The  velocity  vector  is  computed  for  the  combined  velocity 
of  the  earth  and  satellite  and  Herrick's  "f  and  g"  solution  is  used  to 
calculate  the  individual  velocities  each  time  Process  B  accepts  new  data 
from  Process  A. 

We  note  several  points  concerning  the  effects  of  aberration  on  the 
star  tracker.  The  speed  of  the  earth  in  its  orbit  is  30  km/s  and  the 
maximum  speed  of  an  earth  orbiting  satellite  is  8  km/s  relative  to  the 
earth.  Therefore,  the  maximum  shift  in  a  star's  direction  is  about  26 
arcseconds.  This  maximum  occurs  for  stars  90°  from  the  velocity  vector. 
However,  all  stars  in  this  neighborhood  will  be  shifted  by  nearly  this 
amojnt  and,  thus,  the  distortion  of  the  FOV  will  be  insignificant. 

However,  aberration  will  displace  the  boresight  direction.  To  avoid 
orientation  errors  in  the  combined  FOV(A)  and  FOV(B)  solution  we  must 
correct  the  catalogue  direction  cosines  by  applying  equation  (B.6). 

For  those  stars  in  the  direction  of  the  velocity  vector,  the  shift 
in  direction  will  be  small.  But  since  the  shift  is  always  towards  the 
velocity  vector  the  distortion  is  noticeable  (an  apparent  shrinking  of 
the  FOV).  In  this  case,  the  aberration  should  be  applied  before  the  final 
least-squares  solution  for  the  single  FOV. 
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APPENDIX  C:  CALCULATION  OF  ORIENTATION  ERRORS 


We  have  used  Euler  parameters  or  Rodrigues  parameters  to  describe 
the  orientation  of  the  field  of  view  (FOV)  with  respect  to  an  inertial 
frame.  For  our  simulation  studies  we  can  compare  the  calculated  set 
of  parameters  with  the  true  set.  However,  this  comparison  does  not 
reveal  the  effect  the  errors  have  on  the  FOV  orientation.  Therefore, 
we  define  two  angles,  e  and  u,  which  are,  respectively,  the  angular 
error  in  the  calculated  boresight  direction  and  the  rotation  error  about 
the  true  boresight  between  the  calculated  and  true  FOV. 

For  our  calculation  we  consider  only  FOV(A).  Consider  two  reference 

frames  {a^}  and  {ac}  which  are  the  true  and  calculated  frames  for  FOV(A) 

(refer  to  Figure  C-l).  The  magnitude  of  the  cross  product  between  a^  and 
-  c 

a_  ,  the  boresight  vectors,  yields  the  error  e ,  the  angle  between  boresights 
T 

a  *  a„  =  sin  £  -  e. 

C3  T3 

The  calculations  for  angle  p  is  more  involved.  We  first  define  a  new 
frame  {c}  such  that 


c,  =  (a  x  a  ) / sin  c. 

3  c3  T3 

From  this  we  can  determine  the  rotation  matrix  [ CN  ] .  The  unit  vector, 

a  ,  is  rotated  into  the  { c }  frame: 

C1 

d  £  [CN]a  . 

C1 

This  unit  vector  is  rotated  by  the  angle  e  about  c3  in  the  {c}  frame: 
r  :  [R(e)]d. 
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where  [R(e)]  = 


cos  l  sin  e  0 


-sin  e  cos  e  0 


Finally,  this  unit  vector,  r,  is  rotated  back  into  the  inertial  frame 


aR  =  [CN]Tr 


Combining  these  equations  we  write: 


a  =  [CN]  [R(e) ] [CN]a  . 

1 


To  calculate  the  rotation  about  the  true  boresight  we  compute: 


a  x  a 
R  T, 


=  sin  y  -  u 


These  two  angles,  e  and  p,  give  a  good  pointing  error  measure  for  the 
star  tracker  system. 
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